;***************************************************************
;discription： CASA模型主程序入点
;主要完成：程序界面，模型所需影像加载，参数验证，过程、结果等消息打印
;
;author: wangqiang cdut 2015-9-10
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
PRO PRO_CASA_MAIN
  ;严格编译规则
  COMPILE_OPT IDL2
  ;二次开发模式初始化，恢复envi的sav文件
  ENVI,/RESTORE_BASE_SAVE_FILES
  ;初始化envi，保存日志文件
  ENVI_BATCH_INIT,log_file='pro_npp_casa.log'
  ;********************************开始****************************************

  PRINT ,'THE CASA MODEL FOR SIMULATE NPP IS BEGINNING ... ...'
  ;print the begin time
  starttime=systime(1)
  PRINT , '-----' ,starttime, '-----'

  ;read standard envi img file
  ;ndvi_file='D:\Program Files\Exelis\ENVI51\classic\data\can_tmr.img'
  ndvi_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\NDVI_SICHUAN2009ALL_PRJ.tif'
  pre_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\PRE_SICHUAN2009ALL_PRJ.TIF'
  tmp_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\TMP_SICHUAN2009ALL.TIF'
  
  tmpmax_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\TMP_MAX_SICHUAN20080101.TIF'
  tmpmin_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\TMP_MIN_SICHUAN20080101.TIF'
  tmpmean_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\TMP_MEAN_SICHUAN20080101.TIF'
  
  dem_file='D:\workspace\idlprojects\PRO_CASA_MAIN\data_sichuan\New Folder\DEM_SICHUAN_CLIP_PRJ.TIF'
  
  ;open the image file
  envi_open_file, ndvi_file, r_fid=ndvi_fid ; open the ndvi file
  envi_open_file, pre_file, r_fid=pre_fid ; 降水
  envi_open_file, tmp_file, r_fid=tmp_fid ; 气温
  envi_open_file, tmpmax_file, r_fid=tmpmax_fid ; 气温
  envi_open_file, tmpmin_file, r_fid=tmpmin_fid ; 气温
  envi_open_file, tmpmean_file, r_fid=tmpmean_fid ; 气温
  envi_open_file, dem_file, r_fid=dem_fid ; 高程
  ;query the file info
  ;获取行、列号，波段数，数据类型
  envi_file_query, ndvi_fid, ns=nrows, nl=ncols, nb=nbands, Data_Type=data_type, xstart=xstart, ystart=ystart,dims=dims
  print,'the ndvi file info :',ndvi_fid, nrows, ncols, nbands,dims
  
  ;构建NDVI输入矩阵
  dims = [-1, 0, nrows - 1 , 0,ncols - 1]
  CASE data_type OF
    1:BEGIN 
      ndvi_data = BytArr(nrows,ncols,nbands) ; BYTE  Byte
      pre_data = BytArr(nrows,ncols,nbands) ; BYTE  Byte
      tmp_data = BytArr(nrows,ncols,nbands) ; BYTE  Byte
      ;新建保存NPP计算结果的数组
      npp_data = BytArr(nrows,ncols,nbands) ; BYTE  Byte
      END
    2:BEGIN 
      ndvi_data = IntArr(nrows,ncols,nbands)    ;INT  Integer
      pre_data = IntArr(nrows,ncols,nbands)    ;INT  Integer
      tmp_data = IntArr(nrows,ncols,nbands)    ;INT  Integer
      ;NPP计算结果
      npp_data = IntArr(nrows,ncols,nbands)    ;INT  Integer
      END
    3: BEGIN
      ndvi_data = LonArr(nrows,ncols,nbands)    ;  LONG  Longword integer
      pre_data = LonArr(nrows,ncols,nbands)    ;  LONG  Longword integer
      tmp_data = LonArr(nrows,ncols,nbands)    ;  LONG  Longword integer
      ;NPP计算结果
      npp_data = LonArr(nrows,ncols,nbands)    ;  LONG  Longword integer
      END
    4:BEGIN 
      ndvi_data = FltArr(nrows,ncols,nbands)    ;  FLOAT  Floating point
      pre_data = FltArr(nrows,ncols,nbands)    ;  FLOAT  Floating point
      tmp_data = FltArr(nrows,ncols,nbands)    ;  FLOAT  Floating point
   
      ;NPP计算结果
      npp_data= FltArr(nrows,ncols,nbands)    ;  FLOAT  Floating point
      END
    5:BEGIN
      ndvi_data = DblArr(nrows,ncols,nbands)    ;  DOUBLE  Double-precision floating
      pre_data = DblArr(nrows,ncols,nbands)    ;  DOUBLE  Double-precision floating
      tmp_data = DblArr(nrows,ncols,nbands)    ;  DOUBLE  Double-precision floating
      ;NPP计算结果
      npp_data = DblArr(nrows,ncols,nbands)    ;  DOUBLE  Double-precision floating
      END
    6:  ;  COMPLEX  Complex floating
    7:BEGIN
      ndvi_data = StrArr(nrows,ncols,nbands)    ;  STRING  String
      pre_data = StrArr(nrows,ncols,nbands)    ;  STRING  String
      tmp_data = StrArr(nrows,ncols,nbands)    ;  STRING  String
      ;NPP计算结果
      npp_data = StrArr(nrows,ncols,nbands)    ;  STRING  String
      END
    8:  ;  STRUCT  Structure
    9:  ;  DCOMPLEX  Double-precision complex
    10: ; POINTER  Pointer
    11: ; OBJREF  Object reference
    12:BEGIN 
      ndvi_data = uintarr(nrows,ncols,nbands); UINT  Unsigned Integer
      pre_data = uintarr(nrows,ncols,nbands); UINT  Unsigned Integer
      tmp_data = uintarr(nrows,ncols,nbands); UINT  Unsigned Integer
      ;NPP计算结果
      npp_data = uintarr(nrows,ncols,nbands); UINT  Unsigned Integer
      END
    13: ; ULONG  Unsigned Longword Integer
    14: ; LONG64  64-bit Integer
    15: ; ULONG64  Unsigned 64-bit Integer
  ENDCASE

  
  ;获取文件投影，注意所用文件采用wgs84地理坐标系
  map_info = envi_get_map_info(fid=ndvi_fid)
  help,map_info
  print,map_info
  
  
  ;********************************数据准备工作***************************
  ;1、NDVI,降水，气温     
  nodata=-9999.0
  FOR i=0,nbands-1 DO BEGIN
    ;get the image data in the openning file which band pos=i
    ndvi= envi_get_data(fid=ndvi_fid, dims=dims, pos=i) ;波段数从0开始,提取出训练数据,并且第一层数据为目标数据
    ndvi_data[*,*,i]=ndvi[*,*]

    pre= envi_get_data(fid=pre_fid, dims=dims, pos=i) ;
    pre_data[*,*,i]=pre[*,*]

    tmp= envi_get_data(fid=tmp_fid, dims=dims, pos=i) ;
    tmp_data[*,*,i]=tmp[*,*]
  ENDFOR
  
  
  
  ;2、获取高程信息
  dem=envi_get_data(fid=dem_fid,dims=dims,pos=0)
  help,dem
  
  ;3、日最高气温、最低气温、平均气温
  num_days=31
  tmpmax_data = FltArr(nrows,ncols,num_days)    ;  FLOAT  Floating point
  tmpmin_data = FltArr(nrows,ncols,num_days)    ;  FLOAT  Floating point
  tmpmean_data = FltArr(nrows,ncols,num_days)    ;  FLOAT  Floating point
  
  FOR i=0, num_days-1 DO BEGIN
    ;get the image data in the openning file which band pos=i
    tmpmax= envi_get_data(fid=tmpmax_fid, dims=dims, pos=0) ;波段数从0开始,提取出训练数据,并且第一层数据为目标数据
    tmpmax_data[*,*,i]=tmpmax[*,*]

    tmpmin= envi_get_data(fid=tmpmin_fid, dims=dims, pos=0) ;
    tmpmin_data[*,*,i]=tmpmin[*,*]

    tmpmean= envi_get_data(fid=tmpmean_fid, dims=dims, pos=0) ;
    tmpmean_data[*,*,i]=tmpmean[*,*]
  ENDFOR
  
  ;***********************开始逐月进行参数计算*********************************
  ;目前只准备了2009年的数据，所以没有嵌套年的循环，每年按12月循环计算
  ;目前这部分只是计算了一个月份的
  year=2009
  FOR month=0, 11 DO BEGIN
    ;输出文件名称
    outputname='npp_'+strtrim(string(year),2)+'_'+strtrim( string(month+1),2)
    print,strcompress( outputname )
    
  ENDFOR
  
  ;月NDVI
  ndvi_2009xx=ndvi_data[*,*,0]
  ;月总降水量
  pre_2009xx=pre_data[*,*,0]
  ;月平均气温
  tmp_2009xx=tmp_data[*,*,0]
  
  ;npp结果文件
  npp_2009xx=npp_data[*,*,0]
  
  ;********************1参数apar****************************
  ;太阳总辐射量
  sol = FltArr(nrows,ncols,  1)    ;  FLOAT  Floating point
  jmonth=1
  PRO_APAR_SOL,nrows, ncols, dem_fid, tmpmax_data, tmpmin_data, tmpmean_data, jmonth, sol=sol
  ;print,sol 
  
  ;光合有效辐射的吸收比例
  pfar = FltArr(nrows,ncols,1)    ;  FLOAT  Floating point
  PRO_APAR_FPAR,nrows,ncols,ndvi_2009xx,pfar=pfar
  help, pfar
  ;print, pfar

  ;植被吸收的光合有效辐射
  npp_apar = FltArr(nrows,ncols,1)    ;  FLOAT  Floating point
  PRO_NPP_APAR,nrows,ncols,pfar, sol,apar=npp_apar
  help,npp_apar
  ;print,npp_apar
  
  
  ;*********************2参数e***************************
  ;植物生长的最适宜温度 
  PRO_T_OPTIMAL,nrows,ncols, tmp_data,t_opt=t_opt
  help,t_opt
  ;print,t_opt 
  ;低温胁迫因子
  e_t1 = FltArr(nrows,ncols,1)    ;  FLOAT  Floating point
  PRO_E_T1, nrows, ncols, t_opt, tmp_2009xx, e_t1=e_t1
  help,e_t1
  ;print,e_t1
  ;高温胁迫因子
  e_t2 = FltArr(nrows,ncols,1)    ;  FLOAT  Floating point
  PRO_E_T2, nrows, ncols, t_opt, tmp_2009xx, e_t2=e_t2
  help,e_t2
  ;print,e_t2
  
  ;区域实际蒸散发 
  PRO_W_EET,eet=eet
  help,eet
  print,eet
  ;区域潜在蒸散发
  PRO_W_EPT,ept=ept
  help,ept
  print,ept
  ;水分胁迫因子
  PRO_E_W, eet ,ept ,e_w=e_w
  help,e_w
  print,e_w
  ;植被最大光能利用率
  PRO_E_EMAX, e_max=e_max
  help,e_max
  print,e_max
  
  ;实际光能利用率
  npp_e= FltArr(nrows,ncols,1)    ;  FLOAT  Floating point
  PRO_NPP_E,nrows,ncols, e_t1, e_t2, e_w, e_max, e=npp_e
  help, npp_e
  
  
  ;**********************3开始模型运算npp**************************
  ;初级生产力NPP计算
  PRO_CASA_NPP,nrows, ncols, npp_apar, npp_e, npp=npp_2009xx
  help,npp_2009xx
  ;print,npp
  
  
  
  
  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = '输出CASA模型结果...')
  ;ENVI_WRITE_ENVI_FILE, npp, OUT_NAME = output_file,Map_Info = map_info
  ;二进制方式输出
  OPENW,lun, output_file,/get_lun
  ;写入npp结果数据到文件中
  WRITEU,lun,npp_2009xx
  FREE_LUN,lun
  ;写出头文件信息
  ENVI_SETUP_HEAD,fname=output_file, ns=nrows, nl=ncols, nb=1, interleave=0, data_type=data_type, offset=0,map_info=map_info, /write, /open
 

  ;********************************退出***********************************
  ;print the end time
  endtime=systime(1)
  duration=endtime-starttime
  PRINT ,'-----' ,endtime, '-----'
  PRINT ,'THE CASA MODEL ARE RUNNING ' ,duration,' SECOND TIME'
  PRINT ,'THE CASA MODEL FOR SIMULATE NPP IS ENDED ... ...'
  ;关闭envi二次开发模式
  ENVI_BATCH_EXIT
END

